
;;===========================================================
;; Routine to make a simple mass function from the IRDC distances.

conf = omni_read_conffile('./conffiles/survey_info_bgps2.conf')


;; Read probability vectors
restore,'./local/BGPS_V2_pvec.sav',/ver
s = omni_read_cat(conf.cat)

eind = where(constrain.emaf AND constrain.post, nem)
print,nem

pvec      = pvec[eind]
s         = s[eind]
constrain = constrain[eind]


IF ~TAG_EXIST(s,conf.s[0].val,INDEX=iflux) ||  $
   ~TAG_EXIST(s,conf.s[0].err,INDEX=ieflux) THEN BEGIN
   message,'Error: Tags specified in survey-info config file for integrated '+$
           'flux density aperture are not present in the catalog.  Exiting.',$
           /cont
   RETURN
ENDIF




;; mass(m_sun) = 13.1 * d(kpc)^2 * S(Jy)
m = 13.1 * s.(iflux) * 1.0 * pvec.stat.duse[0]*pvec.stat.duse[0] / 1.d6
print,m4_stat(m)


;; Make logarithmic hostogram
plothist,alog10(m),bin=0.15,/noplot,xarr,yarr

;; Calculate EDGES of the histogram bins (1 more than # of bins)
delx = mean(xarr[1:*] - xarr[0:*])
edges = 10.^([min(xarr)-delx,xarr,max(xarr)+delx,max(xarr)+2.*delx] - delx/2.)

;; Add blank bins at both ends (for plotting), and convert xarr back
;; into MASS (rather than LOG(MASS))
xarr = 10.^[min(xarr)-delx,xarr,max(xarr)+delx]
yarr = [0,yarr,0]

;; Calculate DX in real spacem each a little different
dx = edges[1:*] - edges[0:*]
print,n_elements(xarr),n_elements(dx)
print,m4_stat(dx)

myps,'./masses/mass_function.eps',/cmyk

delta = cgSymbol('Delta')
sun = cgSymbol('sun')
xrange = [1d1,1.5d4]
plot,xarr,yarr/dx,psym=10,/xlog,xr=xrange,/xst,thick=3,$
     xtit='Mass [M'+sun+']',/ylog,yr=[2.d-4,1.d1],/yst,$
     ytit=delta+'N/'+delta+'M [M'+sun+'!u-1!n]',xtickformat='exponent10',$
     ytickformat='exponent10',title='Very Rough Preliminary Mass Function'

print,m4_stat(yarr/dx)


FOR j=0,n_elements(edges)-2 DO BEGIN
   IF edges[j] GE xrange[0] AND edges[j] LE xrange[1] THEN BEGIN
      plots,edges[j],10.^(!y.crange[0])
      plots,edges[j],yarr[j]/dx[j],/cont,thick=3
   ENDIF
ENDFOR


oplot,xarr,yarr/dx,psym=2,thick=2,color=cgColor('Cyan')

mfit = alog10(yarr/dx)

mcut = [250.,550.]
ccolor = ['Red','Dark Green']
slope = fltarr(2)

FOR cc = 0,1 DO BEGIN
   
   
   vline,mcut[cc],/xlog,/ylog,color=cgColor(ccolor[cc]),thick=3
   ii = where(xarr GE mcut[cc] AND finite(mfit), nii)
   
   xfit = alog10(xarr[ii])
   yfit = mfit[ii]
   
   res = linfit2(xfit,yfit)
   
   FOR k=0,nii-1 DO print,xfit[k],yfit[k]
   
   print,res
   slope[cc] = res[1]

   oplot,xarr[ii], 10.^(res[0] + res[1]*alog10(xarr[ii])),$ 
         color=cgColor(ccolor[cc]),thick=3
   
ENDFOR

propto = cgSymbol_extra('propto')
al_legend,/top,/right,box=0,textcolor=ccolor,spacing=1.7,$
          ['dN/dM '+propto+' M!u'+string(slope,format="(F0.2)")+'!n']

myps,/done

END
